import pylab as pl
%matplotlib inline
import sys
folder = '/Users/peppe/work'
if folder not in sys.path:
sys.path.append(folder)
import PS4C
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import ascii
import healpy as hp
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
from PS4C.PS4Cast import *
import astropy
from astropy import units as u, constants as C
dir_ps='/Users/peppe/work/PS4C/'
PICASSO inpainting SPASS maps¶
Giuseppe Puglisi , Xiran Bai , Yuuki Omori Sept. 30 ,2019
We implemented the Python Inpainter for Cosmological and AStrophysical SOurces (PICASSO) which is a package aiming at inpaint any kind of maps encoding either Gaussian (e.g. CMB ) or Non-gaussian signals (e.g. Gal. Foregrounds). We made use of several state of the art techniques already adopted for image reconstruction. Most of'em are based on Convolutional Neural Networks (CNN), one of the most promising is the DCGAN . In this posting gonna show the inpainting performances on SPASS TQU maps with DCGAN trained on CMB maps (notice that this isn't logically correct: the training with non-gaussian signal is still ongoing but we oversees to have the trained weights by the end of the week).
How many sources do we expect to detect in the QU SPASS maps?¶
I run point source forecasts with PS4Cast package (Puglisi + 2018) with the SPASS specs, assuming a detection threshold to be $7 \sigma$ ( this threshold is set by the source detection algorithm ) .
nu=2.3
sens= .023
fwhm= 8.9
fsky= 0.3
SPASS = Experiment(ID='SPASS', sensitivity= sens, frequency=nu , fwhm=fwhm , fsky=fsky,nchannels=1,
units_sensitivity='Jy',units_beam='arcmin')
fc =Forecaster(SPASS, ps4c_dir=dir_ps, sigmadetection=7 )
fc.forecast_pi2scaling(verbose=False )
fc ()
fc.print_info()
Notice that the number quoted above are very close to the one reported in Meyers et al. 2017 ( detection of total flux point sources in SPASS survey , no polarized flux). They 've detected about twice the sources we expect, ~23,000 sources since the detection threshold is lower ( $5\sigma$ ). In polarized flux we expect a definitely smaller number of sources ~ 70, with completeness in polarization at $100 $ mJy.
Performing the Source detection algorithm on the SPASS polarization maps¶
Since we are very interested in polarization, i run the Matched Filter in the polarization amplitude map and if a source is detected with signal-to-noise ratio >7 is included in the detection catalogue.
blobmap= hp.read_map('/Users/peppe/work/heavy_maps/Pfiltered_spasswhole.fits', verbose=False )
spassmap =hp.read_map('/Users/peppe/work/heavy_maps/spass_dr1_1902_healpix_Tb.tqu.fits', field=[0,1,2], verbose=False)
wmask = hp.read_map('/Users/peppe/work/heavy_maps/apodized_spass_footprint.fits.gz', verbose=False)
nside=hp.get_nside(wmask)
obspixs = pl.ma.masked_not_equal(wmask,0. ).mask
snrmap =blobmap /blobmap[obspixs].std()
#print (blobmap[obspixs].std() )
exclude_bright_source = pl.where(snrmap[obspixs]<50)[0]
pl.figure(figsize=(15,15))
snrmap =blobmap *wmask /blobmap[obspixs][exclude_bright_source ].std()
#print (( blobmap[obspixs].std() -blobmap[obspixs][exclude_bright_source ].std() ) /blobmap[obspixs].std() )
beampix = hp.nside2resol(nside,arcmin=True)
conv_factor= PS4C.pointsource_utilities.Krj2brightness(freq=2.3*u.GHz,
theta_beam= 8.9 *u.arcmin,
temp= 1. *u.K )
hp.gnomview( pl.sqrt(spassmap[2]**2 +spassmap[1]**2)
, coord=['G'] ,reso=6
, ysize=300 ,max=0.05,
xsize=700, rot=[0,-60] ,sub=311, title='SPASS Polarizat. Amplitude', notext=False, unit='K') ;
obspixs2 = pl.ma.masked_not_equal(wmask,0. ).mask
rms = (blobmap[obspixs2] *conv_factor ).std().value
hp.gnomview( blobmap *conv_factor
, coord=['G'] ,reso=6
, ysize=300 ,
xsize=700, rot=[0,-60],sub=312, min=-5*rms, max=5*rms, cmap=pl.cm.Greys_r, title='Match Filter map ',
unit='mJy',notext=False)
hp.gnomview( snrmap
, coord=['G'] ,reso=6
, ysize=300 ,
xsize=700, rot=[0,-60],sub=313, min=0, max=5, cmap=pl.cm.Blues, title='SNR map ',notext=False)
Constructing the Catalogue¶
from sklearn.neighbors import DistanceMetric
Haversine_distance = DistanceMetric.get_metric('haversine')
def peakfinder(pix, snrmap):
nside= hp.get_nside(snrmap)
tpeak = snrmap[pix].max()
neigh = hp.get_all_neighbours(lonlat=True,nest=False,nside=nside,theta=pix )
while tpeak < snrmap[neigh].max():
tpeak = snrmap[neigh].max()
maskpix =pl.ma.masked_equal(snrmap[neigh], tpeak ).mask
pix, tpeak= peakfinder( neigh[maskpix], snrmap)
#print "peak found", tpeak, pix
return pix,tpeak
def sizefinder(peak , C, r,delta,snrmap ):
r2=hp.nside2resol(nside )*delta
p1 =hp.query_disc( nside, C, r)
p2 =hp.query_disc(nside,C, r2)
ring =pl.array( list ( set(p2)^ set(p1) ))
if snrmap[ring ].mean()<0 :
#print "size found",pl.rad2deg(r)*60
return r
else:
delta = delta +1.
r= sizefinder(peak,C, r2,delta, snrmap )
return r
class Source :
def __init__(self, lon, lat, size , snr, pixel ):
self.lon =lon
self.lat = lat
coords = SkyCoord( l =self.lon *u.deg , b=self.lat *u.deg, frame='galactic')
self.ra = coords.icrs.ra.deg
self.dec = coords.icrs.dec.deg
self.size=size
self.snr =snr
self.pix=pixel
self.extended=False
class SourceList():
def __init__(self) :
self.sources=[]
def add_source (self, source) :
if len(self.sources) ==0 :
self.sources.append(source)
else:
Lat =pl.array( [ s.lat[0] for s in self.sources ] )
Lon =pl.array( [ s.lon[0] for s in self.sources ] )
X=pl.deg2rad( [ (s.lon[0] ,s.lat[0] ) for s in self.sources ] )
pos = pl.deg2rad( (source.lon, source.lat) ).T
sigma = source.size
dmat = Haversine_distance.pairwise(X, pos)
if dmat.min() !=0 :
self.sources.append(source)
def savelist (self, filename):
strings =['GLON' ,'GLAT' ,'RA' , 'DEC' ,'SIZE_ARCMIN', 'SNR', 'EXTENDED', 'MASK_RADIUS' ]
Lat =pl.array( [ s.lat[0] for s in self.sources ] )
Lon =pl.array( [ s.lon[0] for s in self.sources ] )
Ra =pl.array( [ s.ra[0] for s in self.sources ] )
Dec =pl.array( [ s.dec[0] for s in self.sources ] )
size =pl.array( [ s.size *60 for s in self.sources ] )
SNR =pl.array( [ s.snr for s in self.sources ] )
Ext =pl.array( [ s.extended for s in self.sources ] )
vals = [ Lon,Lat, Ra,Dec, size, SNR, Ext, pl.array(self.maskradii ) ]
ascii.write(table = vals,
names= strings, output= filename , overwrite=True )
SNRthreshold=7
blobs=pl.where(snrmap*wmask >= SNRthreshold )[0]
nside= hp.get_nside(snrmap)
mask=pl.ones_like(snrmap)
r1=hp.nside2resol(nside )
beamfwhm =8.9
SList= SourceList()
pmap =pl.sqrt( spassmap[1]**2 +spassmap[2]**2)
for ipix in ( blobs ):
peakpix, snrpeak = peakfinder(ipix,snrmap)
vec = pl.array( hp.pix2vec(ipix=peakpix,nest=False, nside=nside) )
lon, lat = ( hp.vec2ang( vectors=vec, lonlat=True ) )
size =pl.rad2deg( sizefinder(peakpix, vec, r1, 2 ,snrmap) ) #degree
s=Source(lon,lat,size , snrpeak, peakpix )
if size > beamfwhm/60:
s.extended=True
SList.add_source(s)
SList.maskradii=[]
for s in SList.sources:
vec = pl.array( hp.pix2vec(ipix=s.pix,nest=False, nside=nside) )
hole = hp.query_disc(nside=nside, vec=vec, radius=pl.deg2rad( s.size ))
T_integrated= pmap[hole].sum() *u.K
flux =PS4C.pointsource_utilities.Krj2brightness(freq=2.3*u.GHz,
theta_beam= s.size *u.deg,
temp= T_integrated).to(u.Jy)
f = 1. + pl.log10(flux.value)
if s.extended:
radius = pl.deg2rad(s.size *f )
else:
radius = pl.deg2rad(beamfwhm/60. )
SList.maskradii.append(radius)
disc = hp.query_disc(nside=nside, vec=vec, radius=radius )
mask[disc ] = 0.
print ("Found %d point sources w/ SNR> %d" %( len( SList.sources ), SNRthreshold) )
Ext =pl.array( [ s.extended for s in SList.sources ] )
print ("%d are extended " %( len( Ext[Ext] ) ))
SList.savelist('/Users/peppe/work/inpainting/FG_inpainting/SPASS_pol_pointsources_SNR{}.dat'.format(SNRthreshold))
hp.write_map("/Users/peppe/work/heavy_maps/SPASS_pointsource_mask_matchfilter_SNR{}.fits".format(SNRthreshold), mask, overwrite=True )
As you can notice there are ~500 sources most of them are in the Galactic Plane ( in PS4Cast these aren't included ).Thus, we have to apply a mask which selects the Extra-galactic sources (setting a threshold in Galactic latitude $|b|>20 ^{\circ}$, excluding the extended sources ( i.e. the ones whose size is larger than the beam) .
strip = hp.query_strip(nside, pl.deg2rad(70),pl.deg2rad(110) )
strippedmask=wmask*1.
strippedmask[strip]=0.
ones = np.where( strippedmask>0)[0]
hp.mollview(mask*strippedmask)
hp.graticule()
Above the point source mask excluding the low Gal. latitude sources.
tab = ascii.read('/Users/peppe/work/inpainting/FG_inpainting/SPASS_pol_pointsources_SNR7.dat')
Nentries =tab['RA' ].shape[0]
extended = tab['SIZE_ARCMIN'] >=10.
Selecting unresolved sources at $|b|> 20 ^{\circ}$¶
s =0
maps=pl.sqrt( spassmap[2]**2 + spassmap[1]**2)*( wmask )
masklat= pl.ma.masked_greater(abs( tab['GLAT'] ), 20 )
masksize=pl.ma.masked_less_equal(tab['SIZE_ARCMIN'],11 )
maskand = pl.logical_and ( masklat.mask, masksize.mask )
print(tab['GLAT'][maskand].shape)
for lon ,lat in zip(tab['GLON'][maskand] , tab['GLAT' ][maskand]) :
pl.figure()
hp.gnomview(maps , rot=[lon, lat],
xsize=128, reso=3.5,sub=121, min=0 )
hp.gnomview(( snrmap ), rot=[lon,lat],min=0., max=7, xsize=128,
cmap=pl.cm.Blues, reso=3.5,sub=122)
Notice that we end up considering 45 sources, which is very close to our expectations.
Inpainted Maps¶
The masked area is centered in the detection coordinates and a masking radius of 30 arcminutes (3x SPASS beam).
Now let's see what Picasso can do, let's inpaint'em all!
ra,dec= pl.loadtxt(folder+'/inpainting/FG_inpainting/SPASS_sources_snr7_glon_glat.dat')
for i in range(ra.shape[0]):
pl.figure()
im = pl.load(folder+'/inpainting/FG_inpainting/stacks/SPASS/singlestacks/Q_{:.5f}_{:.5f}_masked.npy'.format(ra[i],dec[i] ))
inp = pl.load(folder+'/inpainting/FG_inpainting/outputs/Contextual-Attention/Q_{:.5f}_{:.5f}.npy'.format(ra[i],dec[i] ))
pl.subplot(121)
pl.title('Masked SPASS')
pl.imshow(im);pl.colorbar(orientation='horizontal')
pl.subplot(122)
pl.title('Inpainted SPASS')
pl.imshow(inp);pl.colorbar(orientation='horizontal')
These are 6x6 deg flat sky Q maps (128 x 128 pixels ) extracted from the SPASS map and centered at the source location (sizes are similar to the maps shown above). As already stated above, this can be definitely improved once we 've the trained weights on Synchrotron simulated maps. Overall looks very promising!